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Thermalization has been shown to occur in a number of closed quantum many-body systems, but 
the description of the actual thermalization dynamics is prohibitively complex. Here, we present 
a model - in one and two dimensions - for which we can analytically show that the evolution into 
thermal equilibrium is governed by a Fokker-Planck equation derived from the underlying quantum 
dynamics. Our approach does not rely on a formal distinction of weakly coupled bath and system 
degrees of freedom. The results show that transitions within narrow energy shells lead to a dynamics 
which is dominated by entropy and establishes detailed balance conditions that determine both the 
eventual equilibrium state and the non-equilibrium relaxation to it. 



How a closed quantum system reaches a steady state in 
which observables assume time-independent expectation 
values is both an intriguing and fundamental question in 
physics. When such a steady state can be described by an 
equilibrium thermodynamic ensemble we speak of "ther- 
malization" . The study of closed quantum systems and 
their relaxation behavior has recently received renewed 
interest [THU HHIS] fuelled by ground breaking progress 
in preparing and manipulating ultracold atomic quan- 
tum gases |16| . Current experiments provide excellent 
thermal insulation from the environment and grant co- 
herence times much longer than typical relaxation times. 
These approximately closed systems therefore constitute 
an ideal test bed for the investigation of relaxation phe- 
nomena and thermalization |17H21) . In spite of this ex- 
perimental progress it is fair to say that there is still a 
lack of understanding as to when and under what cir- 
cumstances quantum systems prepared far from thermal 
equilibrium relax to a steady state and when the acquired 
equilibrium state is compatible with predictions from sta- 
tistical mechanics. 

A currently very successful and insightful approach re- 
lates thermalization to spectral properties of the many- 
body Hamiltonian P2l426j . In thermalizing systems the 
expectation value of an observable O calculated in an 
eigenstate with energy e coincides with the thermal aver- 
age taken in the microcanonical ensemble with the same 
energy, i.e. (O)Mc(e)- This eigenstate thermalization hy- 
pothesis (ETH) [22l [22 which states that any individual 
many-body eigenstate contains a thermal state is numer- 
ically verified for a number of systems [H |31 UHl IE] ■ In 
practice any initial state j^*) is a superposition of en- 
ergy eigenstates |a). When these states are chosen to 
lie around a given energy e shell relaxation to the mi- 
crocanonical state will occur for long times as coherences 
between states of different energy wash out due to oscil- 
lating phase factors with incommensurate frequencies j2J . 
While the ETH is very useful to assess whether a system 
thermalises or not, the calculation of the actual thermal- 
ization dynamics requires the solution of the many-body 
Schrodinger equation. The complexity of this task grows 
exponentially with the number of degrees of freedom and 



even the most powerful numerical methods are restricted 
to small systems, low dimensions and/or short propaga- 
tion times (21j . 

From classical statistical mechanics we know that re- 
laxation is often described by an effective evolution equa- 
tion, such as a Master or Fokker-Planck equation. It 
has the eventual equilibrium distribution as its stationary 
state due to the condition of detailed balance of its tran- 
sition rates |27j . It is an intriguing question whether such 
Fokker-Planck equation can be actually also be found for 
a thermalising and closed many-body quantum system. 
If the answer is affirmative it will mean that expectation 
values of O can be calculated not only in the limit of in- 
finite time, but also for intermediate times, and far from 
equilibrium, from a statistical ensemble. 

Finding a general answer to this question is a 
formidable task. However, in order to show that this 
idea can work in principle we present a non-integrable 
spin model in which an effective Fokker-Planck equation 
can be analytically derived - a simple Ising-like spin- 1/2 
model on a lattice in a transverse magnetic field. Its 
Hamiltonian reads 

L L 

H^Hn+Hv = nY,<^k+Vj2nknk+i, (1) 

k=l k=l 

with the Pauh spin matrix cr^ = (|| + ||) and 
projector Uk = (|t) (tl)fe on the up-state on lattice site 
k. The coupling of the spins to the field is parameterized 
by fl. In contrast to the conventional Ising model the 
spin-spin interaction of strength V is state-dependent, 
i.e. two neighboring spins only interact when both of 
them are in the up-state. We study this model because it 
exhibits generic features, such as non-integrability, even 
in one dimension, and local interactions leading to com- 
plex dynamics. Furthermore, it is amenable to analyt- 
ical treatment, and can be realised in experiment with 
ultracold lattice gases of laser-driven Rydberg atoms 
m H El [Ml 130], polar molecules [3T], or with trapped 
ions [SaES]. 

In the strongly interacting regime, V ^ H. this system 
has been shown to thermalize [21 |S] . Thermalization was 
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FIG. 1. a) Graphical representation of the Hilbert space for a system with L = 8 sites and V — )■ oo. Each node represents 
a classical spin configuration with fixed excitation number n. The number of configurations with exciation number near L/4 
grows exponentially in the system size. Links connecting the nodes represent transitions between spin configurations induced 
by the Hamiltonian. We obtain an effective evolution equation of the quantum state on the graph by a second order expansion 
of the von-Neuman-equation. Here only three fundamental moves can be identified: transmissions, reflections and loops. The 
solid colored circles shown represent beginning and end of such move while the hollow circles represent the intermediate state. In 
the thermodynamic limit loops constitute the dominant contribution, b) Properties of the graph, such as transition coefficients 
between columns r„_>„-i-i can be analytically calculated by exploiting a connection to the partition function of a gas of hard 
core dimers in one dimension. An example for a spin configuration together with the corresponding dimer arrangement is 
shown. Deposition/removal of a dimer effectuates a move to the right/left on the graph, c) The dynamics on the graph can be 
mapped to a one-dimensional Fokker-Planck equation describing diffusion of the probability distribution 7r(?y, t) in the potential 
U{y). The variable y is related to the excitation density n/L by a transformation explained in the text. The shape of the 
potential is determined directly by the branching ratio of the graph a). An out-of-equilibrium distribution (yellow) relaxes into 
an equilibrium state (blue) which is centered at the minimum of the potential. Different colors indicated different times in the 
interval < Qt < 1.6. d) Temporal evolution of excitation number distribution Pn{t) for an initial state with 7 excitations in a 
system of 25 sites. The plot shows data obtained from the full quantum calculation (red) in comparison to the result obtained 
from the Master Equation ([5|. Both are in excellent agreement even for short times. 



demonstrated to become explicit in the distribution func- 
tion pn (t) of the number of excitations (up-spins) which 
was shown to reach an equihbrium value for long times, 
i.e. t ^ Pn{t) provides a meaningful measure of 

the degree of thermalization as it does not only capture 
local properties, e.g. the number of excitations 

L L 

(iV)t^^(nfc), = ^np„(i), (2) 

k=l n=0 

but also encodes spatial correlations through 

L L 
k,m=l n=0 

and higher moments. 

To illuminate the thermalization process we consider 
the limit V ^ oo. Here the space of states with finite en- 
ergy is spanned by all spin configurations in which neigh- 
boring excitations are absent. This space can be repre- 
sented as the graph depicted in Fig. flk (shown for L = 8 



for ease of visualization). The nodes of the graph corre- 
spond to classical spin configurations with fixed number, 
n = 0, . . . ,L/2, of (non-contiguous) excitations. Transi- 
tions between the states are driven by the Hamiltonian 
and shown as links connecting the nodes. For example 
the state lititiiii) with n = 2 excitations is directly 
coupled to the configuration lititiiti) with n = 3 up- 
spins but not with litiititi) since both are not linked 
by just a single spin-flip operation. The quantum evolu- 
tion of the system can be thought of as motion within this 
graph, where the amplitudes of a quantum state |4'(i)) 
are related to the occupation of the nodes. 

The distribution function Pnit) is defined as Pn(t) = 
Tr[P„p(t)], where P„ is a projection operator onto the 
subspace spanned by all states contained in the n-th col- 
umn of the graph, see Fig. [T^. The density matrix for 
the whole closed system, p{t) = \^{t)){'i/{t)\, evolves ac- 
cording to the von-Neumann equation (h = 1) 

dtp{t)^^i[H,p{t)]. (4) 
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Starting from Eq. Q it is possible to derive (see Sup- 
plementary material for details) for our spin problem a 
Master equation for the time evolution of the distribution 
function; 

Pn~l{t)] 

-2^(2 t (T„^„_i + T„^„+l) Pnit) (5) 

The crucial insight that leads to this equation is that the 
complexity in the linkage of the graph, Fig. la, allows for 
a statistical treatment of the quantum dynamics. This 
becomes more accurate the larger the system size and 
for initial configurations that are located near the cen- 
tral region of the network, i.e., that are situated not too 
close to the left or right edge. Note that the latter con- 
dition is almost always satisfied for a randomly selected 
initial configuration. Our approach bears similarities to 
the one presented in Ref. [7]. The difference, however, 
is that we do not rely on a system-bath partitioning of 
the quantum system. The transition coefficients T„_j.„±i 
determine the probabilities for creating and annihilating 
excitations. They can be calculated analytically by ex- 
ploiting the fact that the set of spin configurations that 
form the graph, Figjl^, are the same as those of hard-core 
dimers on a one-dimensional lattice, Fig|l|D. The number 
of possibilities to remove/deposit one such dimer from a 
configuration with n dimers is directly related to the typ- 
ical number of links connecting a node in the n-th column 
of the graph. 

Eq. ([5]) is derived using a second order expansion where 
Hamiltonian Eq. ([I]) induces only three kinds of tran- 
sitions in the graph: loop transitions, reflections, and 
transmissions; see Fig. la. In loop transitions the ini- 
tial and final state is identical, while the intermediate 
state is a configuration which differs by one excitation. 
Reflections are generalizations of loop transitions where 
the initial and final states have the same number n of 
excitations but differ in their specific arrangement. For 
transmissions initial and final value of n differ by two. 
Combinatoric arguments (see Supplementary material) 
imply that the contribution of loop transitions to the 
transition coefficients Tn^n±i is by far the most domi- 
nant one. For instance, reflections between two randomly 
chosen states are highly unlikely since due to the struc- 
ture of the Hamiltonian both states must be identical 
except for the position of exactly one dimer. Explicitly, 
the transition coefficients T!„_>„-|-i are. 



Tn- 



-1 = n, 



{L-2n-l){L-2n) 

yn+l - 77 7^ (,Dj 



{L- 



1) 



With these coefficients the Master equation ^ obeys de- 
tailed balance, i.e., p'^^Tn^n^i = p'^^^Tn+i^n for some 
equilibrium distribution p'^. It can be easily verified that 
this distribution is given by [31] 



L 



1 + y/Ej L-n 



L — n 



(7) 



In the continuum limit, L,n ^ 1 we can introduce the 
excitation density x = n/L and the Master Equation ([5| 



becomes a Fokker-Planck equation (FPE) for the distri- 
bution function p{x, t), 



dtp{x, t) 



d^F[x) - -dlD{x) 



p{x,t). (8) 



Notice that this FPE has both drift and diffusion co- 
efficients, F{x) = (1 - 5a; + 5x^)/{l - x) and D{x) ^ 
(1 — 3x + 3a;^)/[(l — x)L], which are dependent on x. 
Note also that while the FPE is local in time it de- 
pends explicitly on t through the overall rate 2n,^t re- 
sulting in a typical time dependence ~ exp [— AfJ^t^] 
with A > 0. This is a direct consequence of the un- 
derlying coherent quantum evolution. The density de- 
pendent diffusion coefficient, D{x), is a manifestation of 
the network structure of the state space in which the 
system moves according to the FPE. Since a; is a contin- 
uous variable, D{x) can be thought of as a metric. All 
one-dimcnsional spaces are flat, so via a suitable coordi- 
nate transformation y = y{x) we can bring the FPE to 
a form dtn{y, t) — —Vt^t[dyF{y) — Ddy]'K{y, t), where the 
new scaled diffusion coefficient is constant, D = 1, and 
P{y) = ~dyU{y) is a gradient force with the potential 
U(y) (see Supplementary material). In the y representa- 
tion relaxation corresponds to the approach to stationar- 
ity in the potential C/(y), see FigjI];. 

A quantitative comparison between the diffusive dy- 
namics predicted by the Master equation ([s]) and the 
exact quantum evolution is given in Fig. [l}i. Here we 
show the temporal evolution of the excitation number 
distribution Pn{t) for a system with L = 25 spins start- 
ing from a pure state with 7 excitations. This state is 
located in the central column of the configuration graph, 
cf. Fig. [1^, where the assumptions that entered in the 
derivation Eq. ^ are best met. That is, the initial con- 
ditions have to correspond to highly connected nodes of 
the graph, so that subsequent dynamics lead to a rapid 
spreading within the network of states. If the initial state 
is near the edge of the graph, for example a state with no 
excitations, its low connectivity implies that the mixing 
is less rapid and coherent oscillations persist for longer 
times [2 . Nevertheless, the vast majority of micro-states, 
and therefore most randomly chosen initial states, are lo- 
cated in the highly connected region of the graph and 
thus evolve according to the Master equation. Indeed, 
Fig. [TJi shows excellent agreement between the Master 
equation and the full quantum evolution for the relax- 
ation of the distribution p„ (i) from an initial state with 
fixed n to the stationary state. Both the thermal dis- 
tribution which is established for il.t > 1.5 and also the 
short time dynamics is captured very accurately by the 
Master equation ([5|. For small values of L there are 
deviations from the thermal-like evolution. These fiuctu- 
ations are finite size effects which we expect to vanish for 
large enough systems: Fig. [2^ shows that as anticipated 
temporal fluctuations become less pronounced with in- 
creasing L. 

The obvious question is now how general these results 
are, for example, do they depend on dimensionality. To 
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FIG. 2. a) Finite size efTects. Shown in black is the evolution 
of the excitation fraction (x) — {N)/L for different system 
sizes, L = 15, 20, 25. Each panel shows data of the quan- 
tum calculation for 5 initial states, randomly chosen from the 
manifold of 3, 5 and 7 excitations. The red curve is the re- 
sult obtained from the Fokker-Planck equation. The temporal 
fluctuations around the Fokker-Planck result decrease with 
increasing system size, b) The results also hold for higher di- 
mensions, here two dimensions. The evolution can also here 
be mapped on a one-dimensional Fokker-Planck equation in 
excitation number space. The corresponding transition coef- 
ficients can be calculated from the partition function of a gas 
of hard squares, c) Equilibrium distribution of the excita- 
tion number (black) for a two-dimensional 6x6 lattice with 
periodic boundary conditions. The agreement with the nu- 
merically exact quantum result (red), here shown at Qt — 3.0 
for an initial state with 8 excitations, is excellent. Note that 
the colormap is the same as in Fig. [TJi. 



address this we have also studied a two-dimensional ver- 
sion of the model ([I]) . Like in the one-dimensional case 
the dynamics can be analysed with the help of a graph. 
Configurations representing the nodes are given by all 
possible arrangements of hard squares on a square lat- 
tice, see Fig. [2]d. While the transition coefficients can 
no longer be calculated analytically, the numerically ob- 
tained Master equation for a 6 x 6 lattice with periodic 
boundary conditions shows a similar correspondence be- 
tween the fully quantum evolution and the approximate 
Master equation as in the one-dimensional case, as illus- 
trated in Fig. [2]:. 



In summary, we have shown for a spin model in one and 
two dimensions that the relaxation of a closed quantum 
many-body system to a thermal steady state is well ap- 
proximated by a Master or Fokker-Planck equation. This 
is derived analytically from the full coherent dynamics, 
for the probability distribution of the thermalizing ob- 
servables. The arguments that lead to such an evolution 
equation are the same that justify the eventual thermal 
steady state: the complex connectivity of the state space 
and the corresponding mixing effectuated by the Hamil- 
tonian. This brings up the question of how the properties 
of the configuration network are related to the statisti- 
cal complexity of the eigenstates that is central for the 
ETH. Clearly, both are connected as they stem from one 
and the same closed system Hamiltonian. Our model is 
physically relevant and yet generic - it is non-integrable, 
has local interactions and no disorder. Yet, an analytical 
analysis of the complex many-body dynamics is possible 
which is unusual for non-integrable systems. It represents 
a concrete example for how a thermalization dynamics 
can emerge from a purely quantum mechanical evolution 
without an explicit system-bath partitioning and/or pos- 
tulated randomness. 
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I. SUPPLEMENTARY MATERIAL 

For the interested reader we present here some tech- 
nical details on the derivation of the Master equation as 
weU as on the transformation of the Fokker-Planck equa- 
tion to a form with constant diffusion term. 



Basis states. In order to reflect the structure of the 
configuration network depicted in Fig. la in our calcu- 
lations, we denote our micro-states (vertices of the net- 
work) in the from |nC„). Here n denotes the number of 
up-spins and C„ their geometrical arrangement (configu- 
ration). The micro-states form a complete, orthonormal 



basis, 



\nC„){nCn\ = 1 



(nC„|m/C„,.) ^ 6 {nCn,mlCm) ■ 



(9) 



Sometimes it will also be convenient to remind us about 
the configuration energy E = {nCn\Hv\nCn) of a micro- 
state. In this case will denote it as \nCn',E). 

Sketch of the derivation of the classical Master 
equation. Using the integral form of the von-Neumann 
equation the time evolution of the expectation value {0)t 
of an observable is given by 

dt{d)t^~ J^dsTr [6 [H lit), [Hi is), Pis)]]}, (10) 

where Hi it) — JJ^ Hq\] is the interaction picture Hamil- 
tonian obtained with the unitary transformation U = 
exp[— itiJy]. In writing eq.(10l we have assumed that 
the initial state is such that the density matrix at f = 
commutes with O. This is satisfied, if we assume that 
our initial state is a micro-state and our observable is 
the number of up-spins with fixed configuration energy 
E, 



O = 



'\nCn\E) inCn\E\. 



(11) 



The prime indicates that the summations has to be taken 
only over configurations with configuration energy E. Af- 
ter inserting eq. (Ill into eq. (10), expanding the double 



commutator and using the cyclic property of the trace to 
interchange the density matrix to the rightmost side of 
each term, we perform three major steps to arrive at the 
final Master equation. 

(i) Since our main interest is in the description of the 
time evolution of initial states that lie deep within our 
complex configuration network (i.e. that are connected 
to a huge number of other states), we assume that the 
system quickly loses its memory about its past. Our first 
approximation is, therefore, to replace pis) — )■ pit) in 
eq. ( |10[ ), which is effectively a second order expansion 
of the von-Neumann equation reminiscent of the Redfied 
equation [1]. 

(ii) We only keep the statistically dominant loop tran- 
sitions in our description (see next section for the justifi- 
cation of this step and ref. |2] for a numerical example). 
This amounts to only keeping the diagonal elements of 
the squared Hamiltonian Hi. After these replacements 
we perform the time integrals that yield functions of 
the form r(t) = sin [iE' - E)t]/ iE' - E). Since we re- 
strict ourselves to the dynamics within one energy shell 
iE' — E) this function simplifies to r(t) t. 

(iii) Using the fact that the Hamiltonian only couples 
states that differ in excitation number by one, we can per- 
form the summation over the excitation number index. 
In our final approximation we replace the remaining ma- 
trix elements of iJ| by times the mean connectivity 
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of the n-manifolds (see below for details) and arrive at 
the Master equation presented in our main manuscript. 

Relative importance of the three transition types. 

In the derivation of the Master equation we only take 
loop-transitions into account (c.f. Fig. la). The neglect 
of the other transition types for n 1 is justified as fol- 
lows: If we randomly pick a state with n excitations the 
number of loop transitions is simply given by the num- 
ber of direct connections of a chosen mirco-state to states 
with one more or less up-spins. This number is on the 
order of magnitude of the system size L. This has to be 
contrasted with the situation faced when regarding reflec- 
tions. If one randomly picks an initial and a distinct final 
micro-state within a manifold having 1 < n < L /2 exci- 
tations, these states will most probably not be connected 
at all by a reflective transition, since the spin-flip term 
of the Hamiltonian can only couple two states that differ 
by the position of one excitation. The probability that 
two states are linked by a reflection can be assessed by 
dividing the number of configuration pairs that exactly 
differ by the position of one up-spin (which is on the or- 
der of L) by the total number of configurations Vn with 
n-excitations (that for n ^ 1 grows exponentially with 
system size). It is vanishingly small for \ < n < L/2. 
By a similar argument one can also see that transmis- 
sions are statistically insignificant as compared to loop 
transitions. 

Calculation of graph connectivity. In order to obtain 
the mean connectivity between manifolds with n — 1 and 
n up-spins we start by determining the total number of 
allowed transitions from states with n to states with n—1 
excitations, 



(12) 



where cr^. is the spin lowering operator at site k and \n) = 
\nCn)- The calculation is most easily done by using 
a spin-coherent state |^) introduced in ref. [5], 



k 

= \Q)-m+em---- 



(13) 



with mfc — 1 — rifc. This state is a weighted sum over all 
allowed micro-states in the E = shell. It is normalized 
to the partition function S(z) of a hard-dimer gas with 
fugacity z = 

Interestingly, expectation values of off-diagonal oper- 
ators taken with the state |^) can be expressed by ex- 
pectation values of diagonal operators. In particular we 
have {^lEk'^klO = -r'(eiEfc"felO with 



1 



1 - 



1 



S(z) 



(14) 



Using these relations we can define a generating function 
A(z) for the coefficients c„_>.„_i 



A(z)=^(e|n,|0=E^"c«-"-i 



(15) 



SO that 



^n-^n— 1 



1 9" 
n! dz" 



K{z) 



z=0 



n i^-f> 



(16) 

The typical number of connections (i.e. possible tran- 
sitions) T„_j.„_i that a randomly ehosen state from the 
manifold with n excitations has to the n—1 -manifold can 



then be obtained by dividing ( 16 ) by the total number 



Vji of states with n excitations generated by ^{z) 

2n-l 



1 9" 



-S(z) 



z=0 



L 

77, ' 



n (L-^y 



(17) 



]=n+l 



Similarly all other coefficients of the Master equation can 
be obtained. 

Fokker-Planck equation and transformation to 
constant diffusion term. Our aim is to transform the 
Fokker Planck equation such, that the diffusive term be- 
comes constant. In general such a transformation reads 

dy 

p{x)dx — 7r(t/)dy = 7r(i/)— dx = jTT{y)dx (18) 
dx 



with the Jacobian J — dy/dx. From (18) we get 



p{x) = J 7r(y) and ^ 



(19) 



In order to obtain constant diffusion we have to require 
J^D = D^const dy{J^D)^Q (20) 

The transformed Fokker-Planck equation then reads 

b . 



9t7r(y,t) 



with 



F = J 



dyF{y) 



Ay,t) (21) 



(22) 



Finally, if we set Z) = 1 it follows from ( 20 ) that 
dy 



J = 



dx 



D{x) 



(23) 



so that the transformation x — > y{x) can in principle be 
explicitly calculated. 



The solution of eq. ( 23 ) can be expressed analytically 



in terms of various elliptic integrals over inverse trigono- 
metric and hyperbolic functions. It is, however, not very 
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illuminating. Luckily, in the range x = ... 1/2 the so- 
lution can excellently be fitted by a quadratic function 
and thus the inverse transformation also be determined 

y{x) = oix + a2x'^ 

<y) = ^ ai ) (24) 



with ai = 0.7074 and 02 = 0.4169. The effective force 
F{y) can now explicitly be calculated. 
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